'''
Purpose:     To compute the Augmented Dickey Fuller test for residuals
             from a conintegrating regression, supporting deterministic 
             polynomial trends.  

Created on Aug 16, 2010
@author: Peter Harrington
peter.b.harrington@gmail.com
'''
from numpy import *
from lag import *
from detrend import *       #there has to be a better way to do this
from tdiff import *         #perhaps after adding pyconometrics 
from trimr import *         #to my path, this won't be needed
from rztcrit import *

def cadf(inMatX, inMatY, p, nlags):
    if (p < -1):
        print "Error: p cannot be less than -1"
    numObs = inMatX.shape[0]
    if ((numObs - 2*nlags + 1) < 1):
        print "Error nlags is too large in cadf, negative degrees of freedom"
        
    inMatX = detrend(inMatX,p)
    inMatY = detrend(inMatY,p)
    b = ((inMatX.transpose()*inMatX).I)*(inMatX.transpose()*inMatY)
    r = inMatY - inMatX*b
    dep = tdiff(r,1)
    dep = trimr(dep,1,0)
    z = trimr(lag(r,1),1,0)
    k = 1
    while (k <= nlags):
        z = concatenate((z,lag(dep,k)),1)
        k += 1
    z = trimr(z,nlags,0)
    dep = trimr(dep,nlags,0)
    
    dtemp = detrend(z,0)
    beta = linalg.solve(dtemp.transpose()*dtemp, dtemp.transpose()*detrend(dep,0))    #beta = a\b
    print "beta is :",beta
    res = detrend(dep,0) - detrend(z,0)*beta
    so = (res.transpose()*res)/(dep.shape[0]-z.shape[1])
    print "the shape of so is: ",so.shape
    print "the shape of z is: ",z.shape
    var_cov = so[0,0] * ((z.transpose()*z).I)
    
    results={}          #use dictonary to return bundled results
    results['alpha'] = beta[0,0]
    results['adf'] = beta[0,0]/sqrt(var_cov[0,0])
    results['crit'] = rztcrit(numObs,inMatX.shape[1],p)
    results['nlag'] = nlags
    results['nvar'] = inMatX.shape[1]   #number of columns in inMatX
    return results